Random Matrix Model for Superconductors in a Magnetic Field 
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We introduce a random matrix ensemble for bulk type-II superconductors in the mixed state 
and determine the single-particle excitation spectrum using random matrix theory. The results 
are compared with planar tunnel junction experiments in PbBi/Ge thin films. More low energy 
states appear than in the Abrikosov-Gor'kov-Maki or Ginzburg-Landau descriptions, consistent 
with observations. 
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Random matrices have been used to understand the 
distribution of level spacings and widths in nuclei E10 , 
complex atoms |2j, small metallic particles and 
quantum systems whose classical limits are chaotic fl. 
Random matrix models have also appeared in the con- 
text of solving certain SU(N)-invariant field theories in 
the large N limit BJtJ and discretizing two dimensional 
quantum gravity S. In this paper we consider a different 
context: using a random matrix model to solve the elec- 
tronic structure problem posed by the BCS description 
of a superconductor in a magnetic field. 

The BCS description for an ideal electronic system, 
and the more general description for a realistic system, 
can be formulated in a particular basis in terms of a large 
matrix in which the matrix elements have rapidly varying 
phases and a smooth magnitude distribution. This prob- 
lem, however, is too complicated to solve exactly. For an 
ideal system, the matrix is too large, and for a realistic 
system, the matrix elements are unknown. There are two 
motivations for using a random matrix model: 1) In the 
limit of small level spacings compared to energy scales 
of interest, spectra are often insensitive to the details of 
matrix elements with uncorrelated phases and the same 
average magnitude. The robustness of the Wigner semi- 
circle distribution is an example. 2) A simple set of inte- 
gral equations may be derived for the spectrum. These 
equations may be solved numerically and compared di- 
rectly with experiments. 

The microscopic model for a superconductor in a mag- 
netic field was developed as a generalization of the zero- 
field BCS H theory by Gor'kov @, using a Green's 
function description, and by de Gennes using a wave- 
function description. The variational Hamiltonian is 
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and the order parameter may be taken to be a constant, 
$(r) = A . Eq. ([!]) then separates into 2x2 matri- 
ces, yielding the BCS spectrum E^ = \J e\ + Aq, where 
e k = h 2 k 2 /2m- E F . 

In magnetic fields large compared to the field of first 
penetration, H c i, the field inside the superconductor is 
nearly uniform. The cigenstates of 7io, for an ideal sys- 
tem, are Landau levels. The momentum with respect 
to the vortex lattice, which distinguishes states within a 
Landau level, is conserved, and Eq. (|l]) separates into 
2N x 2A*" matrices where N ~ ui D /u> c , the cutoff energy 
for the pairing interaction over the Landau level spac- 
ing. This formulation has been used to calculate certain 
electronic properties of an ideal system [|l 2|Jl 3| ] . 

A more general method |l4| , which does not depend on 
the exact eigenstates being Landau levels, is to consider, 
as in the Anderson description of dirty superconductors 
in no magnetic field |jl5| , arbitrary eigenstates ip a of the 
bare Hamiltonian: 
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New quasiparticle operators may be defined with respect 
to this basis: cft aa = Jdr ^ a (r)cJ CT . The order parameter 
may be written 
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where x(r) is normalized so that J |x( r )| 2 dr = O, and 
4>, the spatial average, is real and positive. We can then 
define a pairing matrix 



Aa/3 



so that the variational two-body Hamiltonian is 
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where ^\ = [ cL c r , ] , c\. a are the electron creation op- 
erators, and 7in is the bare Hamiltonian 7in = (l/2m) , , Tl . , c ,. n , Tr + 

' 2 v ' ' where W is the vector of quasiparticle operators w' = 



(iV — eA/c) — E F . The order parameter is 
<P(r) = -V n(c r] c rl ) 
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] and £ is the diagonal matrix of 
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for a system in a volume Q and with a local two-body 
interaction of strength Vq > 0. In the absence of a mag- 
netic field, plane waves diagonalizc the bare Hamiltonian 



bare eigenvalues \e a \ < Hlu d . The cutoff energy for the 
pairing interaction, Tiuj d , is taken to be the Debye energy 
in conventional superconductors. 

In zero magnetic field, time-reversal invariance in an 
isotropic superconductor (\ — 1) insures that the pairing 
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matrix connects only time reversed states: A a p oc 8 a ,p- 
In the case of the ideal system, near H C 2, it has been 
shown that A a p oc exp(— \e a — eg\ 2 /W 2 ) where Wq ~ A , 
the zero-field gap ]l^ , |l3| |. In other words, the effect of 
the magnetic field is to "fuzz" out the energy range over 
which states are paired, by an amount of order the zero- 
field gap. 

The above pairing matrix description |l2]-|l4|] motivates 
the choice of the following model. We consider an ensem- 
ble of 2N x 2N Hcrmitian matrices 
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where Eq is a diagonal matrix of N uniformly distributed 
eigenvalues \si\ < 1, and the pairing matrix A is 



Aij = 7N h{£l 
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where cy are selected from a complex Gaussian random 



distribution (c*jC kl ) 



SikSji. The cutoff function h is 
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The matrix Eq models the spectrum of the normal metal 
in the range E F — Tiuj d < E < E F + huj D ] all ener- 
gies have been normalized in units of the Debye energy. 
The bare energy level spacing 5 = hui D /N determines 
the size of the matrix. The level spacing is assumed to 
be small, so that the N — > oo limit applies. The cut- 
off function h is related to the Fourier transform of the 
time-dependence of the pair correlation function. The 
Lorentzian form used in Eq. (|^) is expected on general 
grounds in a diffusive system (Ref. |Dj ], chap. 8). The 
spectrum is, however, not very sensitive to this choice 
of cutoff function: a Gaussian cutoff yields results which 
differ by < 10%. 

The two parameters of the model are <f> and W. <p 
is the spatial average of the superconducting order pa- 
rameter. It decreases from the gap value in zero field, 
4>(H = 0) = Ao, to zero at the phase transition, 
<j){Hc2) = 0. W is the scale for the relative energy differ- 
ence between electrons being paired. Time-reversal sym- 
metry gives W(H = 0) = 0; the scale increases steadily 
with magnetic field, to W(H C 2) = Wq. Solving the full 
self-consistent equations for a given material would in 
principle determine the best values for (f> and W as a 
function of applied field and properties of that material. 
This would also make the model much more complex. We 
therefore, instead, find solutions for arbitrary (f> and W, 
and consider these as free parameters to be determined 
from the data. 

The spectrum of H may be determined using meth- 
ods similar to those discussed in the context of solving 
large N 1+1 dimensional QCD M or 0+1 dimensional 



matrix- valued <p 4 theory [Q. These techniques were fur- 
ther developed in the context of a large N Anderson dis- 
order model (many orbitals per lattice site) by Wegner 
Jl6| and in analyzing the eigenvalue correlators of various 
matrix models by Brezin and Zee [ p"7| . 

In general, the connection between random matrix 
models and 0+1 dimensional field theory follows from a 
Lagrangian of the form L = Lq + L' , where the bare La- 
grangian contains an N- vector of fermions ip and an TVxiV 
scalar matrix M, Lq = itp^tp + Tr M 2 , and the interaction 
is L' = N~ 1 / 2 i/j^Mtp. In the large N limit, the fermion 
propagator for this theory is the ensemble-averaged 
Green's function for M: J dte iEt (V>lW^ (0)) = ((E - 
M) T. ). This propagator may be evaluated using stan- 
dard Feynman diagram techniques. The important result 
is that in the large N limit only the planar (non-crossing; 
generalized rainbow) diagrams survive |]6|,[7|, ^6|Jl7| l . 

Some additional methods are necessary for the ensem- 
ble defined by Eqs. (0)-@ because of the Pauli matrix 
structure and the cutoff function in A. The Green's func- 
tion for H, G(E) = l/(E — H), may be written in the 
Gor'kov notation, 
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Bold-faced quantities are 2x2 matrices in electron-hole 
space. Dyson's equation is 
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where £ is the self-energy and Go is the bare Green's 
function (<fi — 0). In the large N limit, only the non- 
crossing diagrams survive, which implies 
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where t\ is the Pauli matrix and we introduce the func- 
tion 
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Note that S depends on two energies, E and e. It is 
not a Green's function, it is just a construct useful for 
solving this particular problem. In the non-banded, uni- 
form A case we would have h = 1 and S = (1/N) TrG, 
independent of e. 

Inserting the result for S into Dyson's equation, per- 
forming the sum in Eq. (fl2"|), and replacing the discrete 
bare energies Si by the continuous variable e, yields 
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This is an implicit integral equation for S. Writing 



Si S 3 
S4 S2 



(14) 



and using that Gq 1 = E — r^e, we have that the off- 
diagonal elements of the matrix in the denominator of 
the integral equation are S3 and S4. Hence, a solution 
exists with S3 = S4 = 0. For a given theory, the spectrum 
in the large N limit is unique, so we expect this to be the 
only solution. Eq. (O) then reduces to 
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For a given cj> and W, Eq. ( |15[ ) may be solved numerically 
to obtain Si(e) and S 2 (e) at any energy E. The single- 
particle excitation spectrum then follows from 



Tr G(E) 



dv Si(v,E) 
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and the usual relation p(E) — (— 1/tt) ImTr G{E + ie). 

In the weak coupling limit, the bare BCS gap Ao -C 
huj D , so both </> and W in Eq. (|l5|) are much less than 1. 
In that case, the limits of integration may be extended to 
00 in both directions, W may be scaled out, and solutions 
depend only on the parameter <j)/W . 

In Fig. 1 we compare the spectrum obtained this way 
with tunneling spectra from planar tunnel junctions in 
BiPb/Ge thin films |lq |. These spectra where taken at 
temperatures low enough (T = 360 mK) that thermal 
smearing affects the lineshape by less than a few per- 
cent. This thermal correction is included with the stan- 
dard expression |l^] for the tunneling current dl/dV — 
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FIG. 1. Normalized junction conductance for PbBi/Ge 
with T c = 5.52 K at 360 mK for several magnetic fields [18]. 
Solid lines: random matrix model (4>/W = .3, .8, 1.6). 
Dashed line: Abrikosov-Gor'kov spectrum = .5). 



-(Go/po) JdEp(E)df(E + eV)/8V, where the Fermi 
factor is f(x) = l/(e" x + 1). 

The numerical agreement between the random ma- 
trix model and the tunneling data is fairly close, a lit- 
tle worse at lower fields. Also shown in Fig. 1 is a fit 
with an Abrikosov-Gor'kov spectrum to the 3 T data, in 
which the gap and pair-breaking strength are taken as 
free parameters. A pair-breaking strength j20| C = -5 
was used, which correctly reproduces the peak; a pair- 
breaking strength 1 gives more low energy states but 
fits poorly at all energies. 

One interesting feature of the random matrix model, 
as seen in Fig. 1, is the presence of more low energy states 
at low fields than is conventionally assumed. The conven- 
tional view, however, has been based on several indirect 
arguments rather than a direct solution of a microscopic 
model. 

First, the Abrikosov-Gor'kov solution for the Green's 
functions of a superconductor in the presence of a dilute 
gas of magnetic impurities pl[ | was applied to type-II 
superconductors by Maki, with impurity concentration 
mapping to magnetic field strength pc[ |. This mapping 
assumes, in addition to the dirty limit and local elec- 
trodynamics, translationally invariant Green's functions 
(e.g., spatially uniform order parameter), which can not 
apply to a bulk type-II superconductor at any magnetic 
field above H c %. Further, it requires that the impurity av- 
eraging technique, valid for a dilute gas of weak, uncorre- 
cted, scatterers, apply to the case of a nearly uniformly 
penetrating magnetic field distribution, which seems un- 
likely. 

Second, a vortex solution in the Ginzburg-Landau 
model has a small "core" of size £, the coherence length, 
compared to A, the penetration depth. Although the 
Ginzburg-Landau model does not contain electrons, the 
assumption is that the small core size over which the 
order parameter magnitude is reduced acts as a box in 
which electronic states are confined, leading to a small 
"normal fraction" of bound states. However, the order 
parameter does not act on electrons like a potential in a 
single-particle Schroedinger equation. The order parame- 
ter enters the Bogoliubov-de Gennes equations analogous 
to a spatially varying mass in a Dirac equation. Squaring 
the equation shows that currents also act as a confining 
potential. A more direct analogy might be made with 
a tornado: it is not just objects in the eye of the storm 
which are trapped and move with the tornado, but also 
objects in the circulating currents around the eye of the 
storm. 

Bound states in the currents could explain why both 
the data and the random matrix model indicate the pres- 
ence of more low energy states than expected from a 
Ginzburg-Landau picture. Solutions for the electronic 
structure around an isolated vortex would help address 
this question; initial numerical work shows clear devia- 
tion from Ginzburg-Landau behavior E2j. Note that if 
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the bound states can be shown to extend out to a pene- 
tration depth, then the overlap with neighboring vortices, 
for fields above H cl , would cause a significant bandwidth. 
This would further increase the expected low energy den- 
sity of states. 

Experimentally, a surplus of low energy states is seen in 
the planar tunnel junction measurements ]l8[ ] as well as in 
local STM measurements in the vortex state [E3| and near 



SN junctions 24 1. Deviations from the linear in magnetic 
field heat capacity of the Ginzburg-Landau picture have 
been seen clearly in high-T c superconductors ]25| ] , as well 
as conventional superconductors p6[ . More systematic 
studies of the spectrum would help clarify the issue. 

Another interesting feature of Eq. ( |l5| ) is that there is 
a simple solution in the limit W <C (p. Physically, this 
corresponds to the situation where the energy scale for 
time-reversal symmetry breaking is small, yet many bare 
eigenstates are mixed chaotically. This might apply near 
SN junctions, in strongly anisotropic superconductors, or 
at magnetic fields much smaller than H C 2- (Note that it is 
necessary that many bare eigenstates be mixed, even if W 
is small, so that the large N limit for the matrix ensemble 
applies.) In this small-IF limit, the cutoff function h(x) 
becomes a delta-function, and Eq. ( |l5| ) is solved by: 



(E- 



l-(E 2 - e 2 Y 
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where the solution for a 2 is given by o-i (E, e) — a± (E, — e) 
and all energies scale with <p: E = E/2<fi, e = e/2<fi, 
S{ = (2<j))<7i. The resulting spectrum vanishes linearly 
near zero energy. This offers a possible alternative to the 
d-wave pairing explanation for the cuprate superconduc- 
tors, for which many observations have indicated close to 
a linearly vanishing spectrum j27]]. 

We also note that mesoscopic superconducting dots 
have recently been fabricated |2q| . It is an interesting 
possibility that the matrix ensemble discussed here might 
describe the eigenvalue fluctuations of such systems, just 
as the standard GOE ensemble describes the eigenvalue 
fluctuations of small metallic particles. 
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